#directory, packages and functions 
#set working directory 
setwd("~/Documents/RSfPETCog")


#list packages
pags<-list("openxlsx", "reshape","tidyverse","ggpubr","rstatix", "plyr", "corrplot", "ggplot2","reshape2","lme4","gghalves")
# packages for: loading xlsx files
#long/wide tranform


#load packages
lapply(pags, require, character.only=TRUE)
## Loading required package: openxlsx
## Loading required package: reshape
## Loading required package: tidyverse
## Warning: package 'ggplot2' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand()    masks reshape::expand()
## ✖ dplyr::filter()    masks stats::filter()
## ✖ dplyr::lag()       masks stats::lag()
## ✖ dplyr::rename()    masks reshape::rename()
## ✖ lubridate::stamp() masks reshape::stamp()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: ggpubr
## 
## Loading required package: rstatix
## 
## 
## Attaching package: 'rstatix'
## 
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## 
## Loading required package: plyr
## 
## ------------------------------------------------------------------------------
## 
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## 
## ------------------------------------------------------------------------------
## 
## 
## Attaching package: 'plyr'
## 
## 
## The following objects are masked from 'package:rstatix':
## 
##     desc, mutate
## 
## 
## The following object is masked from 'package:ggpubr':
## 
##     mutate
## 
## 
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## 
## 
## The following object is masked from 'package:purrr':
## 
##     compact
## 
## 
## The following objects are masked from 'package:reshape':
## 
##     rename, round_any
## 
## 
## Loading required package: corrplot
## 
## corrplot 0.92 loaded
## 
## Loading required package: reshape2
## 
## 
## Attaching package: 'reshape2'
## 
## 
## The following object is masked from 'package:tidyr':
## 
##     smiths
## 
## 
## The following objects are masked from 'package:reshape':
## 
##     colsplit, melt, recast
## 
## 
## Loading required package: lme4
## Warning: package 'lme4' was built under R version 4.3.3
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## The following object is masked from 'package:reshape':
## 
##     expand
## Warning in check_dep_version(): ABI version mismatch: 
## lme4 was built with Matrix ABI version 1
## Current Matrix ABI version is 0
## Please re-install lme4 from source or restore original 'Matrix' package
## Loading required package: gghalves
## [[1]]
## [1] TRUE
## 
## [[2]]
## [1] TRUE
## 
## [[3]]
## [1] TRUE
## 
## [[4]]
## [1] TRUE
## 
## [[5]]
## [1] TRUE
## 
## [[6]]
## [1] TRUE
## 
## [[7]]
## [1] TRUE
## 
## [[8]]
## [1] TRUE
## 
## [[9]]
## [1] TRUE
## 
## [[10]]
## [1] TRUE
## 
## [[11]]
## [1] TRUE
##Define functions

#function for returning individual cor.test pvalue per cor in data matrix compatible with corrplot
rcor.pvalue<-function (mat, ...) {
  mat <- data.matrix(mat)
  p <- ncol(mat)
  index <- t(combn(p, 2))
  nindex <- nrow(index)
  p_mat<- matrix(nrow = p, ncol = p)
  
  for (i in 1:nindex) {
    p_mat[index[i,1], index[i,2]]<-cor.test(mat[, index[i,1]], mat[, index[i,2]], ...)$p.value
  }
  
  pMat <- p_mat
  colnames(pMat)<-colnames(mat)
  rownames(pMat)<-colnames(mat)
  diag(pMat)<-0
  
  pMat
}

#define stats parameters to be calculated by aggregate 
myStats<-function( x){
  c(Mean = mean ( x),
    SD = sd( x))
}

#read data
CogNET<-read.xlsx(xlsxFile = "RSNNetworkROIsextracted.xlsx", sheet = "Tabelle1",
                           colNames = FALSE)

colnames(CogNET)<-c("MPFC","LPL","LPR","PCC","SML","SMR","SMS","ACC","AIL","AIR","PFCL","PFCR","SMGL","SMGR","LPFCL","PPCL","LPFCR","PPCR","CA","CP")
#read global mean values from file
my_globals <- read.delim("globals_3mm.txt",header=FALSE, sep=",")

#Intensity normalization, divide uptake in each ROI by global per timeframe
CogNETp<-CogNET/my_globals[1:810,]

#export seed values from MPFC after intensity-norm 
#export<-CogNETp[,c(1:20)]
#write.table(export$MPFC, "Cog_Net_gm.txt", col.names = FALSE,row.names=FALSE)

#define subject variable subject id 
subject_id<-c(rep("RS_fPET001",30), rep("RS_fPET002",30), rep("RS_fPET003",30),rep("RS_fPET004",30), rep("RS_fPET005",30),rep("RS_fPET006",30),rep("RS_fPET007",30),rep("RS_fPET008",30),rep("RS_fPET010",30),
rep("RS_fPET011",30),rep("RS_fPET012",30),rep("RS_fPET014",30),rep("RS_fPET015",30),rep("RS_fPET016",30),
rep("RS_fPET017",30),
rep("RS_fPET018",30),rep("RS_fPET019",30),rep("RS_fPET020",30), rep("RS_fPET021",30),rep("RS_fPET022",30),rep("RS_fPET023",30),rep("RS_fPET024",30),rep("RS_fPET025",30),
              rep("RS_fPET026",30),rep("RS_fPET027",30),rep("RS_fPET028",30),rep("RS_fPET029",30))


#add record id and convert record id into factor
CogNETp$record_id<-subject_id
CogNETp$record_id<-factor(CogNETp$record_id)

#add diagnose to_id gm removed data set 
CogNETp$diagnose[CogNETp$record_id=="RS_fPET001"|CogNETp$record_id=="RS_fPET002"|CogNETp$record_id=="RS_fPET003"|CogNETp$record_id=="RS_fPET005"|CogNETp$record_id=="RS_fPET007"|CogNETp$record_id=="RS_fPET008"|CogNETp$record_id=="RS_fPET010"|CogNETp$record_id=="RS_fPET014"|CogNETp$record_id=="RS_fPET016"|CogNETp$record_id=="RS_fPET022"|CogNETp$record_id=="RS_fPET025"|CogNETp$record_id=="RS_fPET026"|CogNETp$record_id=="RS_fPET028"|CogNETp$record_id=="RS_fPET029"]= "PD"

CogNETp$diagnose[CogNETp$record_id=="RS_fPET004"|CogNETp$record_id=="RS_fPET006"|CogNETp$record_id=="RS_fPET011"|CogNETp$record_id=="RS_fPET012"|CogNETp$record_id=="RS_fPET015"|CogNETp$record_id=="RS_fPET017"|CogNETp$record_id=="RS_fPET018"|CogNETp$record_id=="RS_fPET019"|CogNETp$record_id=="RS_fPET020"|CogNETp$record_id=="RS_fPET021"|CogNETp$record_id=="RS_fPET023"|CogNETp$record_id=="RS_fPET024"|CogNETp$record_id=="RS_fPET027"]= "HC"

#convert diagnose into factor
CogNETp$diagnose<-factor(CogNETp$diagnose)
#show record ids
CogNETp$record_id
##   [1] RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001
##   [7] RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001
##  [13] RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001
##  [19] RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001
##  [25] RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001
##  [31] RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002
##  [37] RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002
##  [43] RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002
##  [49] RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002
##  [55] RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002
##  [61] RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003
##  [67] RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003
##  [73] RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003
##  [79] RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003
##  [85] RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003
##  [91] RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004
##  [97] RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004
## [103] RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004
## [109] RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004
## [115] RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004
## [121] RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005
## [127] RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005
## [133] RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005
## [139] RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005
## [145] RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005
## [151] RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006
## [157] RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006
## [163] RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006
## [169] RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006
## [175] RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006
## [181] RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007
## [187] RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007
## [193] RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007
## [199] RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007
## [205] RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007
## [211] RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008
## [217] RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008
## [223] RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008
## [229] RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008
## [235] RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008
## [241] RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010
## [247] RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010
## [253] RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010
## [259] RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010
## [265] RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010
## [271] RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011
## [277] RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011
## [283] RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011
## [289] RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011
## [295] RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011
## [301] RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012
## [307] RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012
## [313] RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012
## [319] RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012
## [325] RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012
## [331] RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014
## [337] RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014
## [343] RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014
## [349] RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014
## [355] RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014
## [361] RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015
## [367] RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015
## [373] RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015
## [379] RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015
## [385] RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015
## [391] RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016
## [397] RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016
## [403] RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016
## [409] RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016
## [415] RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016
## [421] RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017
## [427] RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017
## [433] RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017
## [439] RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017
## [445] RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017
## [451] RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018
## [457] RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018
## [463] RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018
## [469] RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018
## [475] RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018
## [481] RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019
## [487] RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019
## [493] RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019
## [499] RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019
## [505] RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019
## [511] RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020
## [517] RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020
## [523] RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020
## [529] RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020
## [535] RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020
## [541] RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021
## [547] RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021
## [553] RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021
## [559] RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021
## [565] RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021
## [571] RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022
## [577] RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022
## [583] RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022
## [589] RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022
## [595] RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022
## [601] RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023
## [607] RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023
## [613] RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023
## [619] RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023
## [625] RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023
## [631] RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024
## [637] RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024
## [643] RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024
## [649] RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024
## [655] RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024
## [661] RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025
## [667] RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025
## [673] RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025
## [679] RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025
## [685] RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025
## [691] RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026
## [697] RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026
## [703] RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026
## [709] RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026
## [715] RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026
## [721] RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027
## [727] RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027
## [733] RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027
## [739] RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027
## [745] RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027
## [751] RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028
## [757] RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028
## [763] RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028
## [769] RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028
## [775] RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028
## [781] RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029
## [787] RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029
## [793] RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029
## [799] RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029
## [805] RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029
## 27 Levels: RS_fPET001 RS_fPET002 RS_fPET003 RS_fPET004 ... RS_fPET029
#add timepoints
time<-c(rep(61:90, 27)) #repeat time points 61- 90 16 x
CogNETp$time<-time
#insert MCI based on ID according to RSfPET demographics script
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET003"]="PD-MCI"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET005"]="PD-MCI"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET010"]="PD-MCI"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET014"]="PD-MCI"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET025"]="PD-MCI"

CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET001"]="PD-NC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET002"]="PD-NC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET007"]="PD-NC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET008"]="PD-NC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET016"]="PD-NC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET022"]="PD-NC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET026"]="PD-NC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET028"]="PD-NC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET029"]="PD-NC"

#only non-MCI controls
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET004"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET011"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET012"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET015"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET018"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET019"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET020"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET021"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET023"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET024"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET027"]="HC"

CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET017"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET006"]="HC"
#HC RSfPET017 and RSfPET06 have MCI 

CogNETp$MCI_NC_HC<-factor(CogNETp$MCI_NC_HC)


for (i in 1:ncol(CogNETp[,1:20])) {
  MeanTime <- aggregate(CogNETp[,i] ~ time + MCI_NC_HC, data = CogNETp, FUN = mean)
  MeanTime_df <- cbind(MeanTime[-ncol(MeanTime)], MeanTime[[ncol(MeanTime)]])
  colnames(MeanTime_df)<-c("time","diagnose",colnames(CogNETp)[i])
  for (j in 3:ncol(MeanTime_df)){
    plot <- ggplot(CogNETp, aes(x = time, y = CogNETp[,i])) +
     geom_line(aes(color = MCI_NC_HC, group = record_id), linetype = 3,, linewidth = 1) +
    geom_line(data = MeanTime_df, aes(x = time, y = MeanTime_df[,j], color=diagnose), linewidth = 2) +
      scale_color_manual(values=c("#440154","#BFD200","#1F78B4"))+
      theme_bw()+
      labs(title = "", x = "Time series (Frames)",
           y = "[18F]-FDG Uptake")+
      theme(legend.title = element_blank(),legend.text = element_text(size=20),
            plot.margin=margin(50,20,20,20))+
      theme(axis.text=element_text(size=25),
        axis.title=element_text(size=30,face="bold", margin=margin(t=20)))+
    show(plot)
    ggsave(plot,file=paste0("CogNETp",width=12, height=9,units="in", dpi=400,limitsize=FALSE,colnames(CogNETp)[i],".png"))}}
## function (x, y, ...) 
## UseMethod("plot")
## <bytecode: 0x123535270>
## <environment: namespace:base>
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

## Saving 7 x 5 in image

#calculate variation coefficient 

#calculate variation in subjects
Subjects<-c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029")

#define region names
region<-c("MPFC","LPL","LPR","PCC","SML","SMR","SMS","ACC","AIL","AIR","PFCL","PFCR","SMGL","SMGR","LPFCL","PPCL","LPFCR","PPCR","CA","CP")

#define matrix for variation coefficient 
Vk<-matrix(ncol=27, nrow=20)

#loop over subjects to fill matrix with variataion coefficient per subject 

for(i in 1:length(Subjects)){

Vk[1,i]<- sd(subset(CogNETp,record_id==Subjects[i])$MPFC)/mean(subset(CogNETp,record_id==Subjects[i])$MPFC)*100
Vk[2,i]<- sd(subset(CogNETp,record_id==Subjects[i])$LPL)/mean(subset(CogNETp,record_id==Subjects[i])$LPL)*100
Vk[3,i]<- sd(subset(CogNETp,record_id==Subjects[i])$LPR)/mean(subset(CogNETp,record_id==Subjects[i])$LPR)*100
Vk[4,i]<- sd(subset(CogNETp,record_id==Subjects[i])$PCC)/mean(subset(CogNETp,record_id==Subjects[i])$PCC)*100
Vk[5,i]<- sd(subset(CogNETp,record_id==Subjects[i])$SML)/mean(subset(CogNETp,record_id==Subjects[i])$SML)*100
Vk[6,i]<- sd(subset(CogNETp,record_id==Subjects[i])$SMR)/mean(subset(CogNETp,record_id==Subjects[i])$SMR)*100
Vk[7,i]<- sd(subset(CogNETp,record_id==Subjects[i])$SMS)/mean(subset(CogNETp,record_id==Subjects[i])$SMS)*100
Vk[8,i]<- sd(subset(CogNETp,record_id==Subjects[i])$ACC)/mean(subset(CogNETp,record_id==Subjects[i])$ACC)*100
Vk[9,i]<- sd(subset(CogNETp,record_id==Subjects[i])$AIL)/mean(subset(CogNETp,record_id==Subjects[i])$AIL)*100
Vk[10,i]<- sd(subset(CogNETp,record_id==Subjects[i])$AIR)/mean(subset(CogNETp,record_id==Subjects[i])$AIR)*100
Vk[11,i]<- sd(subset(CogNETp,record_id==Subjects[i])$PFCL)/mean(subset(CogNETp,record_id==Subjects[i])$PFCL)*100
Vk[12,i]<- sd(subset(CogNETp,record_id==Subjects[i])$PFCR)/mean(subset(CogNETp,record_id==Subjects[i])$PFCR)*100
Vk[13,i]<- sd(subset(CogNETp,record_id==Subjects[i])$SMGL)/mean(subset(CogNETp,record_id==Subjects[i])$SMGL)*100
Vk[14,i]<- sd(subset(CogNETp,record_id==Subjects[i])$SMGR)/mean(subset(CogNETp,record_id==Subjects[i])$SMGR)*100
Vk[15,i]<- sd(subset(CogNETp,record_id==Subjects[i])$LPFCL)/mean(subset(CogNETp,record_id==Subjects[i])$LPFCL)*100
Vk[16,i]<- sd(subset(CogNETp,record_id==Subjects[i])$PPCL)/mean(subset(CogNETp,record_id==Subjects[i])$PPCL)*100
Vk[17,i]<- sd(subset(CogNETp,record_id==Subjects[i])$LPFCR)/mean(subset(CogNETp,record_id==Subjects[i])$LPFCR)*100
Vk[18,i]<- sd(subset(CogNETp,record_id==Subjects[i])$PPCR)/mean(subset(CogNETp,record_id==Subjects[i])$PPCR)*100
Vk[19,i]<- sd(subset(CogNETp,record_id==Subjects[i])$CA)/mean(subset(CogNETp,record_id==Subjects[i])$CA)*100
Vk[20,i]<- sd(subset(CogNETp,record_id==Subjects[i])$CP)/mean(subset(CogNETp,record_id==Subjects[i])$CP)*100
}

#transform into dataframe
Vk_data<-data.frame(
  Vk)
colnames(Vk_data)<-Subjects
tVk_data<-t(Vk_data)
colnames(tVk_data)<-region

#in groups HC vs PD 
ID_HCPD<-c("PD","PD","PD","HC","PD","HC","PD","PD","PD","HC","HC","PD","HC","PD","HC","HC","HC","HC","HC","PD","HC","HC","PD","PD","HC","PD","PD")



tVk_group<-as.data.frame(tVk_data)
tVk_group<-data.frame(lapply(tVk_group,as.numeric))
tVk_group$ID_HCPD<-factor(ID_HCPD)
summary(tVk_group$ID_HCPD)
## HC PD 
## 13 14
variables<- colnames(tVk_group[,1:12])
grp<-levels(tVk_group$ID_HCPD)

#function
perm_test<- function(variable, group_var){
  group_var<-as.factor(group_var)
  grp_levels<-levels(group_var)
  
  obs_diff<-tapply(variable,group_var,mean, na.rm=TRUE)
  obs_diff_A<-obs_diff[grp_levels[2]]-obs_diff[grp_levels[1]]
  
 perm_diffsA<-replicate(10000,{
    permuted_groups<-sample(group_var)
    perm_mean_diff_A<-tapply(variable,permuted_groups, mean, na.rm=TRUE)
    perm_mean_diff_A[grp_levels[2]]-perm_mean_diff_A[grp_levels[1]]
  })
  
  

     print(length(perm_diffsA))
 
   
  #p-values
   p_value_A<-mean(abs(perm_diffsA)>=abs(obs_diff_A))
   
   return(c(p_value_A))
 
}



#save results
results<-data.frame(Variable=character(), p_value=numeric(), stringsAsFactors = FALSE)

#loop over variables
for (var in variables){
  p_values<-perm_test(tVk_group[[var]], tVk_group$ID_HCPD)
  
  results<-rbind(results, data.frame(Variable=var, p_value_A=p_values[1]))
}
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
print(results)
##    Variable p_value_A
## 1      MPFC    0.8334
## 2       LPL    0.4866
## 3       LPR    0.2651
## 4       PCC    0.0550
## 5       SML    0.9459
## 6       SMR    0.2244
## 7       SMS    0.0358
## 8       ACC    0.3419
## 9       AIL    0.9251
## 10      AIR    0.6611
## 11     PFCL    0.8232
## 12     PFCR    0.8606
#for multple boxplots we need additional long format 
ggplot(long, aes(x=region, y=value, fill=ID_HCPD)) + 
  geom_half_violin(side="r",trim=FALSE, alpha=0.5, show.legend=FALSE)+
  geom_half_boxplot(side="l")+
  facet_wrap(~region, scale="free")+
 # coord_cartesian(ylim=c(min(long$value), max(long$value)+1))+
   scale_fill_manual(labels=c("HC","PD"),values=c("#1F78B4",  "#BFD200"))+
  theme_minimal()+ 
  labs(title = "", x = "Region",
           y = "VC")+
  theme(text = element_text(size = 15),
         title = element_text(size = 20),
        axis.text.x=element_blank(),
        axis.text.y = element_text(size=9),
        legend.title=element_blank())
## Error in eval(expr, envir, enclos): object 'long' not found
 ggsave("Vk_atlas_boxplots.png")
## Saving 7 x 5 in image
## Error in `geom_line()`:
## ! Problem while computing aesthetics.
## ℹ Error occurred in the 1st layer.
## Caused by error in `[.data.frame`:
## ! undefined columns selected
#with Cog Z from demographic script
load("RSfPET.RData")

#in whole gruop with cog z
plot(tVk_group$SMS,subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$Cog_z)
model<-lm(subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$Cog_z~tVk_group$SMS)
abline(model)

summary(model)
## 
## Call:
## lm(formula = subset(RSfPET, record_id %in% c("RS_fPET001", "RS_fPET002", 
##     "RS_fPET003", "RS_fPET004", "RS_fPET005", "RS_fPET006", "RS_fPET007", 
##     "RS_fPET008", "RS_fPET010", "RS_fPET011", "RS_fPET012", "RS_fPET014", 
##     "RS_fPET015", "RS_fPET016", "RS_fPET017", "RS_fPET018", "RS_fPET019", 
##     "RS_fPET020", "RS_fPET021", "RS_fPET022", "RS_fPET023", "RS_fPET024", 
##     "RS_fPET025", "RS_fPET026", "RS_fPET027", "RS_fPET028", "RS_fPET029"))$Cog_z ~ 
##     tVk_group$SMS)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.77102 -0.35105  0.00295  0.48536  1.92220 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)     1.1496     0.6271   1.833   0.0787 .
## tVk_group$SMS  -0.7714     0.3381  -2.281   0.0313 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8706 on 25 degrees of freedom
## Multiple R-squared:  0.1723, Adjusted R-squared:  0.1392 
## F-statistic: 5.205 on 1 and 25 DF,  p-value: 0.03131
#in whole group with moca
plot(tVk_group$SMS,subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$MoCa_Sum)
model<-lm(subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$MoCa_Sum~tVk_group$SMS)
abline(model)

summary(model)
## 
## Call:
## lm(formula = subset(RSfPET, record_id %in% c("RS_fPET001", "RS_fPET002", 
##     "RS_fPET003", "RS_fPET004", "RS_fPET005", "RS_fPET006", "RS_fPET007", 
##     "RS_fPET008", "RS_fPET010", "RS_fPET011", "RS_fPET012", "RS_fPET014", 
##     "RS_fPET015", "RS_fPET016", "RS_fPET017", "RS_fPET018", "RS_fPET019", 
##     "RS_fPET020", "RS_fPET021", "RS_fPET022", "RS_fPET023", "RS_fPET024", 
##     "RS_fPET025", "RS_fPET026", "RS_fPET027", "RS_fPET028", "RS_fPET029"))$MoCa_Sum ~ 
##     tVk_group$SMS)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.1125  -1.1603  -0.0165   2.0909   4.8614 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     32.912      2.442  13.475 5.76e-13 ***
## tVk_group$SMS   -3.370      1.317  -2.559   0.0169 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.391 on 25 degrees of freedom
## Multiple R-squared:  0.2075, Adjusted R-squared:  0.1759 
## F-statistic: 6.548 on 1 and 25 DF,  p-value: 0.01694
#whole group Cogz scores (demographics file)
ggplot(data=tVk_group,aes(x=SMS, y=subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$Cog_z, color=ID_HCPD, shape=RSfPET$MCI_NC_HC)) +
  geom_point(size=3)+
  
geom_smooth(data=tVk_group,aes(group=1),method='lm', formula= y~x,show.legend = FALSE, color="black", se=FALSE)+
  
    scale_color_manual(labels=c("HC","PD"),values=c("#1F78B4",  "#BFD200"))+
  theme_minimal()+
   theme(legend.position = "none", legend.title=element_blank())+
  labs(title = "", x = "SMS VC",
           y = "Cognition z-score")+
   theme(text = element_text(size = 30),
        axis.text.x = element_text(size=20),
        axis.text.y = element_text(size=20))+
  scale_shape_manual(values=c("HC"=16, "MCI"=17, "PD-NC"=16,"PD-MCI"=17))+
  
geom_text(x=1.5,y=-1.0, label="r = -0.42, p = 0.031", color="black", size=7)
## Warning: The following aesthetics were dropped during statistical transformation: shape.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

 ggsave("Vk_SMS_CogZ.png")
## Saving 7 x 5 in image
## Warning: The following aesthetics were dropped during statistical transformation: shape.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
cor.test(subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$Cog_z,tVk_group$SMS, method="pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  subset(RSfPET, record_id %in% c("RS_fPET001", "RS_fPET002", "RS_fPET003", "RS_fPET004", "RS_fPET005", "RS_fPET006", "RS_fPET007", "RS_fPET008", "RS_fPET010", "RS_fPET011", "RS_fPET012", "RS_fPET014", "RS_fPET015", "RS_fPET016", "RS_fPET017", "RS_fPET018", "RS_fPET019", "RS_fPET020", "RS_fPET021", "RS_fPET022", "RS_fPET023", "RS_fPET024", "RS_fPET025", "RS_fPET026", "RS_fPET027", "RS_fPET028", "RS_fPET029"))$Cog_z and tVk_group$SMS
## t = -2.2814, df = 25, p-value = 0.03131
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.68678067 -0.04166061
## sample estimates:
##        cor 
## -0.4151028
cor.test(subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$MoCa_Sum,tVk_group$SMS, method="spearman")
## Warning in cor.test.default(subset(RSfPET, record_id %in% c("RS_fPET001", :
## Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  subset(RSfPET, record_id %in% c("RS_fPET001", "RS_fPET002", "RS_fPET003", "RS_fPET004", "RS_fPET005", "RS_fPET006", "RS_fPET007", "RS_fPET008", "RS_fPET010", "RS_fPET011", "RS_fPET012", "RS_fPET014", "RS_fPET015", "RS_fPET016", "RS_fPET017", "RS_fPET018", "RS_fPET019", "RS_fPET020", "RS_fPET021", "RS_fPET022", "RS_fPET023", "RS_fPET024", "RS_fPET025", "RS_fPET026", "RS_fPET027", "RS_fPET028", "RS_fPET029"))$MoCa_Sum and tVk_group$SMS
## S = 4576.6, p-value = 0.04033
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## -0.397